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\  ABSTRACT 

A  method  for  solving  the  quasineutral  hybrid  plasma  equations  in 
two  dimensions  is  presented,  using  full  ion  dynamics  and  inertialess 
electrons.  The  method  is  extended  to  allow  plasma-vacuum  interfaces 
of  arbitrary  shape.  The  algorithm  is  applied  to  the  study  of  rotational 
instabilities  in  theta  pinch  Vlasov  equilibria. 
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I.  INTRODUCTION 

Many  macroscopic  problems  in  plasma  physics  are  characterized  by  ion  Larmor 
radii  comparable  to  the  scale  lengths  of  the  system.  For  these  problems,  and  for  prob¬ 
lems  involving  microinstabilities,  a  fluid  description  of  the  ions  is  inadequate,  and  the 
ions  must  instead  be  treated  in  a  fully  kinetic  manner.  Also,  as  plasma  behavior  dom¬ 
inated  by  ion  physics  generally  evolves  on  time  scales  much  longer  than  characteristic 
electron  time  scales,  for  many  problems  it  is  unnecessary  to  follow  the  full  dynamics 
of  the  electron  motion.  In  addition,  when  the  frequencies  of  interest  are  low  com¬ 
pared  to  the  ion  cyclotron  frequency,  a><ucl,  the  effects  of  high  frequency 
phenomena,  such  as  electromagnetic  radiation  and  electron  inertia,  are  generally  negli¬ 
gible.  These  considerations  have  led  to  the  development  of  quasineutral  hybrid 
plasma  simulation  codes.  Such  one-dimensional  hybrid  simulations  have  been  per¬ 
formed  by  Sgro  and  Nielson1  for  them  pinch  implosions  and  by  Byers  et  al.2  for  the 
study  of  microinstabilities.  Two-dimensional  (r-z)  simulations  of  theta  pinch  implo¬ 
sions  have  been  performed  by  Hewett3. 

Our  hybrid  algorithm  treats  ions  as  particles  and  electrons  as  an  inertialess  fluid. 
The  Darwin  version  of  Maxwell’s  equations  G.e.,  neglecting  the  transverse  displace¬ 
ment  current)  is  used.  The  electron  momentum  equation,  with  inertial  terms 
neglected,  is  coupled  with  Maxwell’s  equations  and  the  statement  of  quasineutrality  in 
order  to  determine  the  electric  and  magnetic  fields.  Section  II  of  this  paper  describes 
our  two-dimensional  quasineutral  model  which  has  been  applied  to  ion  layer  kink 
instabilities.  Section  in  discusses  the  extension  of  this  algorithm  to  problems  having 
plasma-vacuum  interfaces.  In  Sec.  IV  we  show  the  application  of  the  method  to  the 
study  of  theta  pinch  rotational  instabilities. 


II.  MODEL 


This  section  describes  the  quasineutral  model  and  the  basic  algorithm  used  in 
simulation.  Ampere’s  law  may  be  decomposed  into  its  longitudinal  (curl-free)  and 
transverse  (divergence-free)  parts: 


Vx3-±=I+J.l5. 

c  c  dt 

0  4*7,  |  i  af, 


(la) 


(lb) 


c  c  dt 

where  the  subscripts  /  and  r,  respectively,  refer  to  the  longitudinal  and  transverse 


parts  of  a  vector  quantity.  We  assume  quasineutrality,  setting  which  implies 

that  V  /— 0.  If  7,  vanishes  at  the  boundaries,  or  if  the  system  has  periodic  boun¬ 
daries,  then  7,-0  throughout  the  system.  For  the  study  of  low  frequency  phenomena. 


the  Darwin  approximation  is  made,  that  is,  the  transverse  displacement  current  is 
neglected.  Ampere’s  law  then  reduces  to  the  simple  form 


Vxg-i2L(7,+7/)  (2) 

c 

where  7t  and  7t  refer  to  the  electron  and  ion  current  densities. 

Electrons  are  treated  as  a  fluid  so  that  their  motion  is  assumed  to  be  described  by 
the  electron  momentum  equation 


nt  m,  -  -  enf  (£+7,  x  £/  c)  -  V  Pt  (3) 

dt 

where  me  is  the  electron  mass,  the  electron  density,  v,  the  electron  drift  velocity, 
and  Pe  the  scalar  electron  pressure.  For  low  frequency  modes  electron  inertia  effects 
are  not  important.  Therefore,  the  left  hand  side  of  Eq.  3  is  set  equal  to  zero.  With 
the  electron  current  expressed  as  7,— envf  and  the  assumption  of  quasineutrality, 
Eqs.  2  and  3  may  be  combined  to  produce  an  expression  for  the  electric  field 
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(Vx5)x5--!-7,x5--Lv(/i(7;)  (4) 

4ir/ij€  n,ec  n,e 

where  Te  is  the  electron  temperature.  Equation  4  determines  the  electric  field  as  a 
function  of  the  magnetic  field,  electron  temperature,  ion  current,  and  ion  density.  Ion 
currents  and  densities  are  determined  by  linear  weighting  (particle-in-cell).  The  mag* 
netic  field  is  advanced  by  Faraday's  law 

d3/df--cVx£  (5) 

and  the  ion  particle  positions  and  velocities  are  determined  by  integrating  the  equa¬ 

tions  of  motion 

*/<*--*- (2+ vx3/c).  (6a) 

m, 

die/  dt—~v  (6b) 

A  two-dimensional  simulation  code  has  been  developed  using  this  model.  The 
equations  are  solved  in  cartesian  coordinates  with  no  variation  in  the  z-direction  (i.e., 
d/dr— 0).  The  magnetic  field  is  B—Bzz  and  the  quantities  Jz,  Ez,  and  vz  are  all  set 
equal  to  zero.  Particle  motion  is  followed  by  a  standard  leap-frog  integrator.  The  time 
advance  of  field  quantities  given  by  Eqs.  4  and  5  is  accomplished  by  a  predictor- 
corrector  algorithm.  This  algorithm,  as  follows,  is  similar  in  form  to  one  used  in  one¬ 
dimensional  computations  by  Byers  et  al2.  If  the  quantities  3f+^,  v"  ,  2", 

and  3"  are  known,  the  magnetic  field  is  advanced  by 

3'T+T-3',-(cAf/2)Vx2\  (7a) 

A  prediction  is  then  made  for  2"+I  and  B”*1  by 

+j_ 

Ig — ?n+2lQinlXPt)  2  (7b) 

3£,l-2  2-(cAr/2)VxZ^ 


(7c) 
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Using  the  predicted  fields,  a  predictor  particle  move  is  performed  to  obtain  «/pr£  and 
7 TjL  aft*1  which  is  predicted  by 

-  B$ -  (cA t/2)  Vxl^1.  (7d) 

Finally,  the  new  electric  and  magnetic  fields  are  obtained  from 

En+l-  ai,nltB,Pt)n+T+  ^EOhntXPt)^  (7e) 

S'*1-®  2  —  (cAf/2)  Vx2"+I.  (70 

The  particle  positions  can  now  be  advanced  to  n+y  using  these  new  field  quantities. 

The  electric  and  magnetic  fields  are  stored  on  interlaced  grids.  Spatial  derivatives  are 
determined  by  four-point  operators,  e.g., 

(dEJdx)  i  i  "  [(£,)  ,+UJ+l+(Ex) h-i ,j~(Ex)  ,j+t-(Ex)  U1  /  2  Ax  (8) 

2  'J  2 

This  algorithm  is  second  order  accurate  in  time  and  space.  The  time  step  is  limited  by 
a  Courant-Friedrichs-Lewy  (CFL)  condition  on  the  Alfven  speed,  Ar<(Ax/v*), 
where  vA = B0c/  (4ir  n,  m,) 1/2.  An  additional  constraint  on  the  time  step  for  stability  is 
oiclAr<2;  however,  in  practice  the  CFL  condition  is  the  more  restrictive  requirement 
The  algorithm  of  Eqs.  7  has  been  successfully  applied  in  a  code  using  doubly  periodic 
boundaries  to  study  kink  instabilities  in  field-reversed  ion  layers;  these  results  are  dis¬ 
cussed  elsewhere4. 

IH.  PLASMA- VACUUM  INTERFACES 

The  preceding  algorithm  has  the  disadvantage  that  it  cannot  treat  properly  low 
density  or  vacuum  regions,  where  n,— 0.  We  desire  to  avoid  the  details  of  sheath 
regions  or  describing  low  density  plasma  regions  with  high  accuracy,  since  the 
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quasineutral  hybrid  model  is  not  well  suited  for  such  problems.  However,  we  would 
like  to  include  the  gross  effects  due  to  the  fields  extending  into  vacuum  regions,  such 
as  wall  stabilization.  Additionally,  in  highly  nonlinear  problems,  density  fluctuations 
may  occur  which  cause  low  density  or  vacuum  (n,— *0)  regions  to  arise  in  a  small 
number  of  isolated  cells.  These  will  cause  local  violations  of  the  CFL  condition  which 
are  sufficient  to  terminate  the  simulation  unless  an  alternate  method  can  be  found  for 
determining  the  field  quantities  in  these  cells.  One  possible  solution  is  the  addition  of 
a  low  density  plasma  throughout  the  vacuum  region  which  maintains  sufficient  plasma 
density  so  that  the  CFL  condition  is  satisfied.  However,  this  method  does  not  pro¬ 
duce  the  instantaneous  signal  propagation  that  should  occur  across  a  vacuum.  Addi¬ 
tionally,  for  a  highly  nonlinear  problem,  large  amplitude  waves  are  likely  to  occur  in 
the  low  density  region  which  can  still  cause  the  CFL  condition  to  be  violated. 

It  is  desirable  to  use  a  method  to  treat  vacuum  regions  that  does  not  involve  the 
monitoring  of  complicated,  moving,  plasma-vacuum  interfaces.  One  such  method  is 
described  by  Hewett3,  where  the  resistivity  is  varied  in  moving  from  the  plasma  to  the 
vacuum  regions  in  a  nonlinear  ADI  solution.  We  have  developed  another  way  to 
include  vacuum  regions  by  modification  of  Eqs.  7b  and  7e  in  the  algorithm  of  Sec.  II. 
Two  constants,  cp  and  c„  are  defined  so  that  in  the  vacuum  c,-0  and  cr-l  while  in 
the  plasma  cp- 1  and  cv*0.  Then  Eq.  7b  is  replaced  by 

+! 

c,2£j  +cvV2f£j  -  Cp  (-?"  +  2£(7„  a„ B,Pe)  *  2 )  (9a) 

and  Eq.  7e  is  replaced  by 

c,IM  +  c vV’Z-'-c,  <9b) 

A  grid  point  is  defined  to  be  plasma  if  a,.  and  vacuum  if  n,<  ,  where  ^  is  a 
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cutoff  density,  as  illustrated  in  Fig.  1.  The  effect  of  the  modification  is  that  Eqs.  9a 
and  9b  solve  the  plasma  equations,  Eqs.  7b  and  7e,  in  the  plasma  and  then  solve 
V22!—0,  the  vacuum  electric  field  solution,  in  the  vacuum.  The  magnetic  field  is  still 
advanced  by  Eqs.  7a,  7c,  7d,  and  7f.  These  equations  provide  values  for  Bz  in  the 
plasma,  in  the  vacuum,  and  at  the  conducting  wall.  The  vacuum  field  is  then 
smoothed  by  solving  V2ffz-0  in  the  vacuum  with  the  known  values  of  Bz  in  the 
plasma  and  along  the  conducting  wall  as  boundary  conditions.  The  left  hand  side  of 
Eq.  9  is  constantly  changing  as  the  plasma  moves;  however,  it  is  solved  easily  by  an 
iterative  matrix  solver.  Note  that  no  monitoring  of  the  plasma-vacuum  interface  is 
required.  The  Courant-Friedrichs-Lewy  condition  is  now  imposed  on  the  Alfven 
speed  at  the  density  cutoff,  ne.  For  the  magnetic  field,  this  algorithm  effectively  solves 
V2/T-0  in  the  vacuum.  This  method  adequately  treats  fluctuations  in  low  density 
plasma  regions  and  is  capable  of  handling  highly  dynamical  problems,  as  well.  How¬ 
ever,  for  the  study  of  slowly  evolving  plasma  behavior,  such  as  theta  pinch  rotational 
instabilities  with  growth  rates  much  less  than  the  ion  cyclotron  frequency,  additional 
care  must  be  taken,  as  follows. 

In  our  two-dimensional  model  the  theta  pinch  equilibrium  is  a  cylindrical  plasma 
with  an  azimuthal  current  density  7 -/,(/■)$.  The  initial  plasma  density  is  finite  at  radii 
r<  rp.  A  vacuum  region  extends  from  r—rp  to  a  conducting  wall  at  r— /■*.  An  axisym- 
metric  theta  pinch  equilibrium  must  satisfy  the  radial  force  balance  equation 


(10) 


where  0/  is  the  mean  ion  rotational  frequency.  Near  rp,  the  first  term  on  the  right 
hand  side  of  Eq.  10,  the  magnetic  force  term,  represents  the  inward  force  due  to  mag¬ 


netic  pressure.  The  other  two  terms,  the  thermal  pressure  term  and  the  centrifugal 
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force  term,  are  both  outward  near  rp.  An  equilibrium  may  easily  be  set  up  which 
obeys  Eq.  10  over  the  bulk  of  the  plasma.  However,  near  the  plasma-vacuum  inter¬ 
face  the  algorithm  of  Eqs.  7  and  9  produces  diffusion.  The  reason  for  this  diffusion  is 
that  in  a  region  with  plasma  density  below  the  cutoff  density,  nc,  the  magnetic  field  is 
constant,  because  this  region  is  considered  to  be  vacuum  and  V2A— 0  is  effectively 
solved  there  (rather  than  4tt7/c).  Since  the  current  is  set  to  zero  where 
nc>  n,>0,  the  magnetic  force  term  in  this  region  is  also  zero.  Consequently,  the 
plasma  near  the  interface  will  diffuse  radially  because  it  only  sees  the  outward  forces 
due  to  the  plasma  pressure  and  the  centrifugal  force  of  the  rotating  ions.  The  error 
here  corresponds  to  the  neglect  of  the  small  amount  of  current  density  (both  ion  and 
electron)  carried  by  the  low  density  plasma,  which  in  the  laboratory  provides  the  7x  B 
force  that  contains  the  plasma  at  its  edge. 

An  alternate  method  for  treating  plasma-vacuum  interfaces  has  been  devised  for 
our  two-dimensional  simulation  code.  The  plasma  is  now  considered  to  be  in  three 
regions,  rather  than  two.  The  regions  are  shown  schematically  in  Fig.  2.  The  plasma 
region  is  still  defined  to  be  that  with  density  above  nc.  The  vacuum  region  is  defined 
to  be  that  with  zero  density  or  n,<  n„  where  nv  is  a  cutoff  value  smaller  than  nr.  The 
transition  region  is  the  region  with  ny<n,<  nc.  The  electric  field  is  solved  as  described 
by  Eq.  9.  The  magnetic  field  in  the  plasma  region  is  determined  by  the  plasma  equa¬ 
tions,  Eqs.  7.  However,  to  maintain  the  current  below  the  cutoff,  n^.,  a  different 
vacuum  solution  is  used.  For  the  vacuum  region  the  plasma  is  assumed  to  be  a  per¬ 
fect  conductor,  hence,  vacuum  magnetic  flux  is  assumed  to  be  conserved.  For  this 
two-dimensional  problem  the  magnetic  field  at  a  given  time  will  be  uniform 
throughout  the  vacuum.  Therefore,  B!*  may  be  determined  by 
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B^(t)~  fBfKt-QbdA/A^U) 

vac 

where  Aw  is  the  area  of  the  vacuum  region.  Once  the  vacuum  and  plasma  solutions 
are  known,  the  two  solutions  are  connected  by  solving  V2B,-0  in  the  transition 
region  with  Dirichlet  boundary  conditions  given  by  the  magnetic  fields  determined  in 
the  vacuum  and  plasma  regions.  The  effect  of  this  solution  on  the  equilibrium  is  that 
the  magnetic  field  now  corresponds  to  that  which  would  be  present  if  the  missing 
current  below  the  cutoff  existed  in  the  transition  region.  If  3/8®— 0,  this  current  in 
the  transition  region  is  distributed  so  that  the  transition  azimuthal  current  density  is 
proportional  to  1/r,  since  V2i?,— 0  is  solved  in  that  region.  The  transition  region  is 
typically  only  one  or  two  grid  cells  in  thickness;  however,  maintaining  a  finite  bBjbr 
there  is  critical,  as  our  results  demonstrate.  Sometimes  small  internal  vacuum  regions 
may  occur  due  to  fluctuations  which  arise  in  the  nonlinear  stage  of  a  rotational  insta¬ 
bility.  It  is  possible  to  search  for  such  small  regions  and  to  handle  them  as  transition 
cells  by  solving  V2Br-0.  The  method  does  not  yet  appear  to  be  easily  applicable  to 
complicated  plasma-vacuum  interfaces  in  which  large  internal  vacuum  regions  exist. 

IV.  APPLICATION  TO  THETA  PINCH  ROTATIONAL  INSTABILITIES 

We  have  studied  rotating  theta  pinches  starting  from  exponential  rigid  rotor 
Vlasov  equilibria5.  With  ft,  defined  as  the  mean  rotational  frequency  of  the  elec¬ 
trons,  the  current  density  as  a  function  of  radius  is  described  by 

Ir2±r2 

— j—  (12) 

and  the  corresponding  magnetic  field  profile  is 
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The  sign  preceding  r}  in  Eqs.  12  and  13  is  positive  for  non-reversed  equilibria  (i.e., 
when  5,(0)>0).  For  0,-0  or  r{  sufficiently  large,  a  negative  sign  preceding  rf 
corresponds  to  a  field-reversed  equilibrium  (i.e.,  B2(0)< 0).  A  parameter,  a,  is 
defined  such  that  a— ft  ,/(fl,-0().  The  case  a-0  corresponds  to  a  stationary  ther¬ 
mal  (Maxwellian)  ion  distribution,  with  all  of  the  current  carried  by  the  electrons. 
This  configuration  should  be  stable,  according  to  finite  Larmor  radius  fluid  theory6  and 
Vlasov-fluid  theory7.  We  now  demonstrate  the  difference  between  the  two-region  and 
three-region  solutions.  Simulations  have  been  performed  with  a  density  cutoff,  nc,  at 
three  per  cent  of  the  peak  density.  The  spatial  grid  is  a  100  by  100  mesh.  50,000  parti¬ 
cles  are  used  to  represent  the  ions  and  the  time  step  is  o*CIAr— 0.1.  Figure  3  shows  the 
evolution  of  the  radial  density  profile  (averaged  over  0)  with  time  for  the  two-region 
solution.  The  plasma  diffusion  at  the  edge  eventually  leads  to  the  diffusion  of  the 
bulk  of  the  plasma,  demonstrating  that  with  the  two-region  solution  it  is  impossible  to 
maintain  an  equilibrium  for  times  long  compared  to  the  ion-cyclotron  period.  In  con¬ 
trast,  the  result  for  the  three-region  solution  is  shown  in  Fig.  4.  The  diffusion  at  the 
edge  has  been  virtually  eliminated  and  the  equilibrium  density  profile  is  maintained.  It 
was  found  possible  to  reduce  the  plasma  expansion  in  the  two-region  solution  by  using 
a  lower  cutoff  density;  however,  since  the  time  step  is  limited  by  the  CFL  condition 
on  the  cutoff  density,  unacceptably  small  time  steps,  from  an  economics  standpoint, 
were  required  to  eliminate  the  expansion  effectively.  It  may  be  possible  to  use  a  more 
implicit  method,  for  which  the  time  step  is  not  restricted  by  the  CFL  condition  at  the 
cutoff  density,  in  order  to  improve  the  results  of  the  two-region  solution. 

We  have  done  simulations  starting  from  non-reversed  rigid  rotor  theta  pinch 
equilibria  with  a>0.  These  simulations  have  (ft,-ft /)/«f/-0.022  and  /Sq— 0.75.  0O  is 


the  "beta-on-axis",  0o™/»/(O)  TI(B^ISir)t  where  n,(0)  is  the  initial  density  at  r— 0  and 
Bo  is  the  magnitude  of  the  external  magnetic  field.  No  instabilities  have  been 
observed  for  non-reversed  equilibria  with  a— 1.0.  However,  for  a— 2.0  an  m— 2  insta¬ 
bility,  where  m  is  the  azimuthal  mode  number,  is  evident.  These  observations  are  in 
agreement  with  the  theoretical  predictions  of  Freidberg  and  Pearlstein6  and  Seyier7. 
Figure  5  shows  the  initial  ion  particle  positions  for  an  equilibrium  with  a— 2.0.  Seyier7 
predicts  the  growth  rate  for  this  equilibrium  to  be  (y/o*c/)— 0.025  and  the  real  fre¬ 
quency  to  be  iurlioci)— 0.033.  The  two-region  solution  is  not  feasible  for  this  problem 
because  the  plasma  expansion  would  obscure  the  instability  on  these  time  scales. 
Therefore,  the  three-region  method  is  used.  In  our  simulation,  at  /-144<u~1,  an  m-2 
instability,  where  m  is  the  azimuthal  mode  number,  has  grown  to  large  amplitude  as 
can  be  seen  from  the  ion  particle  positions  shown  in  Fig.  6. 

The  quantity  8r(/,0)3<r>— <r(f— 0)>  is  stored  as  a  diagnostic  and  decom¬ 
posed  into  its  azimuthal  Fourier  components.  In  Fig.  7  is  plotted  (Sr)2  for  mode  2, 
which  gives  an  approximate  measure  of  instability  growth.  The  broken  line  in  Fig.  7 
corresponds  to  (Sr)2  for  the  growth  rate  of  Seyier7  and  shows  that  the  growth  we 
observe  in  our  simulation  is  comparable.  We  have  estimated  the  real  part  of  the  fre¬ 
quency,  mr,  by  measuring  the  rotation  of  the  elliptically  deformed  plasma  cross- 
section.  For  this  simulation  run  we  find  (wr/«CI)-0.035  ±  0.005.  The  value  deter¬ 
mined  by  Seyier7,  (cu,/a>d)— 0.033,  is  within  the  range  of  error  of  our  simulation  esti¬ 
mate. 

Further  details  of  the  code  are  presented  in  Appendix  A.l. 
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V.  SUMMARY 

A  two-dimensional  predictor-corrector  method  for  quasineutral  plasma  simula¬ 
tion,  similar  to  that  used  in  one-dimensional  computations  by  Byers  et  al.2,  has  been 
developed.  The  method  has  been  extended  to  allow  the  inclusion  of  vacuum  and  low 
density  regions.  A  simple,  two-region,  treatment  of  the  plasma-vacuum  interface  has 
been  found  inadequate  for  the  study  of  instabilities  with  growth  rates  much  smaller 
than  the  ion-cyclotron  frequency.  For  such  problems  we  have  implemented  a  three- 
region  method  of  solution  which  avoids  the  problems  of  diffusion  at  the  interface 
found  in  the  two-region  method  and  avoids  finding  the  details  of  the  sheath.  The 
simulation  model  has  been  successfully  applied  to  the  study  of  ion  layer  kink  instabili¬ 
ties  and  theta  pinch  rotational  instabilities. 


•  •  iK.% 
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FIG.  1.  A  schematic  of  the  grid  used  in  the  two-region  solution.  The  shaded  area 
corresponds  to  the  plasma  region,  n,  >  nc  and  the  unshaded  region  to  the  vacuum, 
n,<  n(.  In  order  to  calculate  the  field  quantities  at  a  given  point,  the  operator  to  be 
applied  is  determined  automatically  by  whether  that  grid  point  is  in  the  vacuum  (v)  or 
in  the  plasma  region  (p). 
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A.1  Description  of  Simulation  Code 


The  algorithm  used  to  advance  the  field  quantities  in  our  two-dimensional  hybrid 
code,  AQUARIUS,  is  similar  to  a  one-dimensional  method  used  by  Byers  et  al.1  Both 
methods  have  similar  stability  characteristics.  The  corrector  iteration  is  necessary  to 
prevent  spurious  odd-even  oscillations  from  producing  undesirable  high  frequency 
noise.  The  improvement  with  the  corrector  iteration  can  be  seen  in  Figs.  1  and  2 
which  show  the  oscillation  of  Alfven  waves.  Figure  1,  the  result  without  a  corrector 
iteration,  clearly  shows  the  odd-even  oscillation.  This  odd-even  oscillation  is  not 
present  in  the  predictor-corrector  results  of  Fig.  2.  These  figures  are  taken  from  our 
results  using  a  one-dimensional  version  of  AQUARIUS.  We  have  found  convergence 
to  be  sufficient  after  one  correction,  so  that  further  iterations  are  unnecessary. 

Unlike  the  algorithms  of  Byers  et  al.1,  we  use  B  directly  from  Faraday’s  law 
rather  than  computing  it  from  the  vector  potential,  A.  This  choice  was  made  because 
in  two  dimensions,  without  axisymmetry,  direct  computation  avoids  the  difficulty  of 
decomposing  the  electric  field  into  its  longitudinal  (curl-free)  and  transverse 
(divergence-free)  parts  for  use  in  Faraday’s  law.  We  have  chosen  to  calculate  ?  and 
3  on  separate  grids  because  it  prevents  the  numerical  diffusion  and  instabilities  often 
associated  with  the  computation  of  both  quantities  on  the  same  grid.  This  choice  is 
similar  to  the  procedure  used  in  other  electromagnetic  codes;  e.g.,  Boris  et  aL2.  Our 
predictor-corrector  method  appears  to  be  simpler  than  the  two-dimensional  nonlinear 
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ADI  method  of  Hewett3.  However,  the  more  implicit  nature  of  Hewett’s  method  may 
allow  a  lower  density  cutoff  for  a  given  time  step. 

Our  simulation  model  consists  of  four  separate  parts: 

1.  RREQUI  is  a  program  which  loads  particles  into  an  exponential  rigid  rotor 
distribution  in  r,  0,  v„  and  v*  The  coordinates  are  then  transformed  to  x, 
y,  vx,  and  vy  and  written  to  a  file,  AQLOAD,  as  input  to  AQUARIUS. 

2.  AQUARIUS  is  the  main  program  which  solves  the  time  dependent  field  and 
particle  equations. 

3.  AQPPMD  is  a  post-processor  which  reads  as  input  a  file  produced  by 
AQUARIUS.  This  input  file  contains  the  average  radii  of  the  plasma  for  64 
different  angles  at  selected  time  intervals.  The  post-processor  then  performs 
a  Fourier  decomposition  in  8  to  determine  history  information  concerning 
individual  modes.  This  history  information  is  written  in  the  form  of  a 
graphics  file. 

4.  AQPPHST  is  a  general  post-processor  for  history  information.  This  post¬ 
processor  reads  a  history  file  from  AQUARIUS  which  contains  an  assort¬ 
ment  of  information  at  selected  time  intervals,  such  as  particle  trajectories, 
kinetic  energy,  and  field  energies.  This  information  is  plotted  and  written  to 
a  graphics  file. 

AQUARIUS  simulations  begin  by  the  reading  of  two  input  files: 

1.  INAQ  provides  approximately  SO  basic  parameters  for  the  simulation,  such 
as  grid  size,  time  step,  and  type  and  frequency  of  diagnostics. 


.4 
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2.  AQLOAD  provides  the  particle  positions,  as  well  as  parameters  to  determine 
the  magnetic  held  initialization  for  rigid  rotor  distributions  as  calculated  by 
RREQUI. 

After  the  input  data  has  been  read  and  the  initial  particle  positions  and  magnetic  field 
quantities  determined,  the  electric  field  is  initialized  by  the  equation 

£-  -i_(Vx2?)xS-  -1-7x5-  —VPr 
Avne  nee  ne 

The  initialization  is  finally  completed  by  the  subroutine  SETV  which  performs  a  back¬ 
ward  integration  of  the  ion  particle  velocities  by  one  half  time  step,  so  that  the  leap¬ 
frog  mover  will  be  properly  centered  with  positions  and  velocities  one  half  time  step 
out  of  phase. 

Once  the  initialization  is  complete,  the  main  time  loop  begins.  This  loop  is 
repeated  until  IT,  the  time  step  number,  equals  NT,  the  designated  number  of  time 
steps  for  the  run.  Each  time  step  begins  with  a  particle  move  using  the  subroutine 
MOVER.  MOVER  performs  a  standard  leap-frog  advance  of  the  particle  positions  and 
velocities.  The  subroutine  MOVER  also  accumulates  ion  densities  and  ion  current 
densities  by  linear  weighting  (particle-in-cell).  A  thorough  discussion  of  particle-in-cell 
techniques  and  leap-frog  particle  movers  can  be  found  in  Ref.  4.  MOVER  is  vector¬ 
ized  for  implementation  on  a  Cray-1  computer.  The  vectorization  has  been  found  to 
produce  a  factor  of  2.7  increase  in  speed  over  the  same  routine  without  vectorization. 
After  the  particle  move  is  completed,  the  subroutine  FIELDS  is  called.  FIELDS  per¬ 
forms  the  field  advance  described  in  the  main  section.  It  utilizes  the  subroutines 
CALCE  and  BADV  to  advance  the  electron  and  magnetic  fields.  When  low  density  or 
vacuum  regions  are  present,  CALCE  and  BADV  call  subroutines  to  perform  matrix 
solutions  to  Eqs.  9a  and  9b.  The  matrix  solutions  use  successive  over-relaxation  with 


points  ordered  so  that  this  method  becomes  fully  vectorizable.  An  iterative  scheme 
was  chosen  because: 

1.  As  the  plasma  moves  the  vacuum  regions  change,  changing  the  matrix 
structure  each  time  step,  so  that  a  direct  method  would  have  to  be  fully 
repeated  each  time  step. 

2.  A  good  guess  is  available  from  the  previous  time  step. 

3.  High  accuracy  for  the  vacuum  fields  is  not  required. 

Successive  over-relaxation  was  chosen  because  it  may  be  vectorized  easily  and  because 
it  requires  no  additional  storage  of  temporary  vectors. 

After  the  field  calculation  has  been  completed  in  the  time  loop,  diagnostic  infor¬ 
mation  is  calculated.  In  the  input  data,  different  'snapshot*  diagnostics  may  be 
requested  at  selected  time  intervals,  such  as: 

1.  Phase  space  plots. 

2.  Field,  current,  and  density  contour  plots. 

3.  Field,  current,  and  density  cross-section  plots. 

4.  Bz(r),  £,(/■),  Ee(r),  and  n,(r)  averaged  over  theta. 

5.  B2(9),  Er(9),  E9(9),  Ju(9\  Ji<9(9),  and  n,(0)  for  a  selected  radius. 

6.  Line  printer  contour  plots. 

Additionally,  history  information  is  written  at  selected  intervals  to  serve  as  input  for 
the  post-processors  AQPPMD  and  AQPPHST.  After  the  last  time  step,  IT  "NT,  the 
output  files  are  closed  and  AQUARIUS  ends.  Because  all  history  information  is  writ¬ 
ten  into  post-processor  files,  it  is  unnecessary  to  complete  a  given  run  in  order  to 


access  the  information. 


The  dominant  part  of  the  time  required  to  run  AQUARIUS  is  used  in  MOVER. 
Each  particle  move  takes  approximately  Sfis  per  time  step  of  Cray- 1  CPU  time.  When 
there  are  no  vacuum  regions,  the  time  taken  by  the  field  solve  is  negligible.  For  a  typ¬ 
ical  theta  pinch  run,  with  large  vacuum  regions,  the  field  solve  generally  takes  about 
0.2  s  for  a  100  by  100  grid;  this  is  about  half  of  the  time  required  to  move  50,000  par- 
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